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ABSTRACT 

We extend our previous studies of shock waves and shock-free solutions in thin 
accretion and winds in pseudo-Newtonian geometry to the case when the flow is 
two-dimensional and around a Kerr black hole. We present equations for fully general 
relativistic viscous transonic flows and classify the parameter space according to 
whether or not shocks form in an inviscid flow. We discuss the behaviors of shear, 
angular momentum distribution, heating and cooling in viscous flows. We obtain a 
very significant result: we find that in weak viscosity limit the presence of the standing 
shock waves is more generic in the sense that flows away from the equatorial plane can 
produce shock waves in a wider range of parameter space. Similar conclusion also holds 
when the angular momentum of the black hole is increased. Generally, our conclusions 
regarding the shape of the shock waves are found to agree with results of the existing 
numerical simulations of the two dimensional accretion in Schwarzschild geometry. In 
a strong viscosity limit, the shocks may be located farther out or disappear completely 
as in the pseudo-Newtonian geometry. 
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1. INTRODUCTION 

The recognition that the accretion process onto compact objects are inherently three- 
dimansional is as old as the subject of accretion itself. The pioneering work of Bondi (1952), 
however, contained the solution for spherically symmetric inflow, thereby reducing the problem in 
one dimension. Subsequently, Lynden-Bell (1978) considered the flow to have a constant angular 
momentum but without any radial velocity. He considered the axisymmetric case, so that the 
problem is two dimensional, but still integrable because of the choice of negligible radial motion. 
This study was later extended to understand thick accretion flow structures in pseudo-Newtonian 
geometry (e.g., Paczyhski & Wiita, 1980) as well as in full general relativity (Chakrabarti 1985, 
hereafter C85 and references therein) using more general power law distribution of angular 
momentum which mimicked a viscous rotating flow. 

Nonetheless, major progress on the theoretical front has been done mostly for accretion flows 
which are thin and confined close to the equatorial plane. The study of optically thick viscous 
Keplerian disks (Shakura & Sunyaev, 1973; Novikov & Thorne, 1973), and viscous transonic 
disks (Paczyhski & Bisnovatyi- Kogan, 1981; Chakrabarti, 1990, hereafter C90; for a general 
discussion, see, Chakrabarti 1996a, 1996b, hereafter C96b, and references therein) are some of 
the examples. Secondly, most of these studies are done using the so-called pseudo-Newtonian 
potential of Paczyhski Sz Wiita (1980) which mimics the external geometry of Schwarzschild black 
hole quite well, although the computational error could be significant very close to the horizon 
since the inner boundary condition is not exactly satisfied, (see, Yang and Kafatos, 1995, for a 
discussion on isothermal disks in Schwarzschild geometry) Of interest to us in the context of the 
present paper is a significant observation from the study of the most general viscous transonic 
flow that in a large range of parameter space, spanned by the two free parameters, namely the 
specific energy and angular momentum, the flow develops standing shock waves close to a black 
hole (C90, C96b). Questions naturally arise: Do these shocks continue to be present when a full 
general relativistic study is made? If so, do the general dependences of shock behavior on flow 
parameters remain similar as in the pseudo-Newtonian treatment? Would these shock waves also 
form away from the equatorial plane when the flow is two dimensional? How would the shock 
locations change as one goes from a non-rotating black hole to a rotating black hole? In the present 
paper, we shall answer these questions. Recent works (Chakrabarti & Wiita, 1992; Chakrabarti 
& Titarchuk, 1995) indicate that shock waves could be a natural sites of the hot radiation in a 
black hole accretion, and that the observed quasi-periodic oscillations could be due to some type 
of resonance oscillation in the post-shock flow (Molteni, Sponholz & Chakrabarti, 1996) or due to 
dynamical oscillation of the sub-Keplerian flow (Ryu, Chakrabarti & Molteni, 1996). These make 
it all the more important to understand the behavior of shock waves in accretion and wind flows 
close to a black hole. Of course, a large region of the parameter space do not form shocks at all. 
We therefore classify the entire parameter space for a conical wedge shaped flow to show which 
regions form shock waves, and which regions do not. Another topic of importance which did not 
fetch adequate attention in the literature so far, (despite being very relevant for the self-consistent 
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treatment of a viscous transonic flow) is: how should one treat the inner edge of the flow near the 
horizon? Traditionally, shear is treated to be zero at the inner edge (at r = r ms , the marginally 
stable orbit for a Shakura-Sunyaev [1973] disk) or on the horizon for pseudo-Newtonian viscous 
transonic flows (C90, C96b). Does shear indeed vanish on the horizon? If not, what effect, if any, 
would it have on the angular momentum distribution, or heating processes? We plan to answer 
these questions as well. The motivation is to build a general framework on which viscous transonic 
flows may be studied in future with full consistency when additional heating and cooling effects 
are introduced. 

Numerical simulations of fully two dimensional accretion around a Schwarzschild black hole 
have already provided partial answers to some of the questions posed above. Hawley, Smarr 
& Wilson (1984, 1985, hereafter HSW84, HSW85; Hawley, 1984, hereafter H84) showed that 
quasi-spherical rotating flows develop traveling accretion shock waves whose shapes generally 
follow the funnel walls around the vertical axis. More recent works of Molteni, Lanzafame 
& Chakrabarti, 1994, hereafter MLC94; Ryu et al, 1995; Molteni, Ryu & Chakrabarti, 1996, 
hereafter MRC96) are examples of numerical simulations which produced standing shock waves in 
thin as well as thick accretion flows. The importance of these latter works is that the codes have 
been tested quite adequately against the corresponding analytical transonic solutions (which were 
non-existent during the simulations of HSW84, HSW85) for thin flows, and therefore the existence 
of standing shocks in thicker flows which they reveal may be trusted. 

In the following sections, we shall try to solve the equations of viscous transonic flows in an 
axisymmetric, thick flow from a purely analytical point of view. To achieve our goal we make 
simplifying approximations. We also use fully general relativistic equations for the flow. Accreting 
matter in this case automatically fulfills the inner boundary condition on the horizon. We obtain 
solutions with and without shocks, both for accretion and winds. We arrive at several important 
conclusions regarding shock properties in a thick flow: (a) the region of flow parameters which 
enables shock formation is much bigger in a rotating black hole. However, the general behavior 
of the parameter space itself is similar to what is obtained in in pseudo-Newtonian geometries 
(Chakrabarti 1989, hereafter C89; C90; Sponholz & Molteni, 1994) and (b) We find that the 
shock waves form away from the equatorial plane more easily because of weaker gravity. Since 
these shocks are centrifugally supported, a knowledge of the angular momentum distribution is 
important. We therefore also explore the nature of the shear stress which governs the transport 
of angular momentum. Viscosity not only transports angular momentum, it heats up the gas. 
Assuming the viscosity is due to ions only we compute the heating rate of the flow both in the 
subsonic and the supersonic branches. The cooling is governed by the density and temperature of 
the flow. We present the variation of the bremsstrahlung cooling rate along the flow as well, so 
that we may find out the relative importance of heating and cooling effects which might affect the 
flow properties. 

The organization of the paper is as follows: In Section 2, we present the basic mathematical 
equations for the study of viscous, transonic flows in a black hole geometry. In Section 3, we 
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present a 'brute force' method to find the general behavior of the global surfaces, such as the funnel 
wall, centrifugal barrier, and the shock surfaces. In Section 4, we present the fully self-consistent 
solution of the steady flow which may or may not include shock waves. We then extend the result 
to a thicker flow. In Section 5, we summarize our work and make concluding remarks. 



2. MODEL EQUATIONS 

In what follows, we choose units of velocity, distance, and time to be c, GM^/c 2 and GM^/c 3 
respectively, where G and are the gravitational constant and the mass of the black hole 
respectively. The matter distribution is assumed to be described by the stress-energy tensor 
which for a perfect fluid takes the form, 

= pU^Uy + P{g (lv + Ufj,U v ), (1) 

which satisfies the equation of motion, 

VT^ = 0, (2) 

Here, p is the isotropic thermal pressure and p = po(l + tt) is the mass density, tt being the internal 
energy. The four velocity components satisfy the normalization relation, 

«X = -1 (3) 

where p = 0,1, 2, 3. We assume the metric around a Kerr black hole in Boyer-Lindquist coordinate 
(e.g., NT73), 

ds 2 = g^dx^dx" = - 1 —-dt 2 + -~(d(/> - uodt) 2 + ^—dr 2 + dz 2 (4) 



Where, 



A = r 4 + r 2 a 2 + 2ra 2 

A = r 2 - 2r + a 2 
_ 2ar 
U ~ ~A~ 

Here, g^ u is the metric coefficient and is the four velocity components: 

1/2 



and 

where, the angular velocity is 



A 



(l-V 2 )(l-m)(g H + lg t4) ) 
U(j> = -lu t 



q = ^_ = 9t<p + htt 
u f g<h<h + 
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the I is the specific angular momentum. The radial velocity V in the rotating frame is 

V = 



{i -my/ 2 

where, 

,_Ur^ )V 2 

u t u l 

Since by definition, V = 1 on the horizon and sound speed is expected to be less than unity even 
for the extreme equation of states, the black hole accretion has to be supersonic (and therefore 
sub-Keplerian) on the horizon. Thus, every flow must deviate from a Keplerian disk close to a 
black hole. The equation for the balance of the radial momentum is obtained from, 



which is, 



where, 



(u„u v + 9llv )Tg = 0, (5) 



dr r A r 6 r b r l p + pdr 



7 2 = [l-TZi(^--)]^ 



$ = u r , 

A 2 

Ar 4 

B = (Qa- l) 2 - n 2 r 3 . 

Here and hereafter we use a comma to denote an ordinary derivative and a semi-colon to denote a 
co-variant derivative. The baryon number conservation equation (continuity equation) is obtained 
from 

= 0, 

which is, 

M = 27rrtf£ = 2itA&po (7) 

where, A is the cross section area of the flow. The equation of the conservation of angular 
momentum is obtained from 

= 



or, 



where, 



Pou^ihu^)^ = (rial);-/ (8) 



is the coefficient of dynamical viscosity and v is the coefficient of kinematic viscosity. When 
rotation is dominant the relevant shear tensor component cr^ is given by (Anderson & Lemos, 
1988; Straumann, 1991), 

** = 2^ (9) 
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so that the angular momentum equation takes the form, 

£-£+ = -_L^^/2 3 A l/2 (1Q) 

Where, 

C = —hu^ 

and h is the specific enthalpy: h = (p + p)/po- Here we have modified earlier works of NT73 and 
Peitz (1994) in that the fluid angular momentum C, rather than the particle angular momentum 
—us is used. That way, for an inviscid flow (rj = 0) one recovers £=constant as in a fluid picture. 
Similarly, the radial velocity term is included (eq. 6) and angular momentum is allowed to be 
non-Keplerian (Eq. 10). is the angular momentum on the horizon since shear (as defined 
by Eq. 9) vanishes there. In presence of significant radial velocity, the shear in eq. (9) is to be 
replaced by it's full expression, 

where P^ v = g^ u + u^u 1 " is the projection tensor and = is the expansion (e.g., Shapiro & 
Teukolsky, 1983). In that case, £ + will no longer be the specific angular momentum of the flow 
at the horizon, but at the place where the residual shear effect vanishes. If, on the other hand, 
viscosity is allowed to vanish sufficient rapidly (as we assume in what follows) on the horizon, then 
it could be the angular momentum on the horizon as well. 

Entropy generation equation is obtained from the first law of thermodynamics along with 
baryon conservation equation: 

= Q + -Q- = [2 W ^]/T - Q-, (11) 

which in expanded form, becomes, 

dh _ 1_^ = 2j/S v _ 
dr po dr 

where Q + and Q~ are the heat generation rate (vertically integrated) and the heat loss rate 
respectively, h is the specific enthalpy: h = (p + p)/po- Here, we ignore the terms due to radiative 
transfer. Using rotational shear as given in eq. (9), the entropy equation takes the form, 

dr po dr r b 

Along with the equations given above, one needs to solve the relativistic shock conditions if 
the flow has a standing shock. A shock surface is a timelike three surface formally defined as, 

C(^) = (14) 
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across which some of the thermodynamic variables suffer a first order discontinuity (Taub 1978; 
Anile & Russo 1986). The mass and the energy- momentum conservation across this discontinuity 
demands that 

° C(pu")_ = ^C(pnn + (15«) 



and 



9 AT*")- = ^C(Tn + (156) 



dx^ dxt 1 

where, the signs — and + denote the quantities before and after the shock respectively. The above 
equations could be written as, 

p_n M n M = p+u^n^ (16a) 

and 

p_n u + (p_ + p^)(u^n il )u u _ = p+n u + (p + + p + )(u^n^)u+ (166) 

Here, n M is the four normal vector component across the shock. 

In the next two sections we shall provide explicit examples of global solutions of the above set 
of equations, with or without shocks. First, we present an approximate method to obtain global 
shock surfaces. 



3. SURFACES OF STRONG SHOCKS IN ACCRETION AND WINDS 

From Eq. (16b) one notices that at a shock surface, the sum of ram pressure (due to radial 
motion) and thermal pressure is continuous across the surface. Presently, we wish to obtain the 
nature of this surface assuming that in the in the strong shock limit one can ignore the pre-shock 
thermal pressure compared to the ram pressure. In this limit, (Novikov & Thorne, 1973), 

7 + 1 

P+ tP- (17a) 

7-1 

and 

v+ ~ ^— (176) 

7 + 1 

where, 7 is the adiabatic index. Strictly speaking above assumption is valid for isothermal flows 
7 = 1 only. However, to get a rough approximation of the global behavior, we assume that this 
is valid for 7 = 4/3 also. For the post-shock quantities, we ignore radial motion to compute 
thermal pressure, although we include the ram pressure obtained using Eq. (17b) along with the 
thermal pressure to obtain the shock surfaces. (We do not include turbulent pressure here since it 
involves an extra parameter, although we believe it's effect is important as concluded by MLC94.) 
Thus, assuming that the post-shock region is predominantly rotating, one obtains from Eq. (2) 
(e.g.,Bardeen, 1973), 

p + p u t ~ {i-my 1 j 
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In this limit the equation could be integrated for any power law angular momentum distribution 

(C85), 

I = (19) 

provided a barotropic equation of state p = p(p) is chosen. Here, and denote two positive 
constants with O.O<J,n0<^O.5 in the disk region. This latter constant actually represents the 
viscosity of the flow: for a small viscosity, angular momentum is nearly constant and ~ 0. 
For high viscosity, the angular momentum is nearly Keplerian, n<f, ~ 0.5. Also, the von-Zeipel 
parameter (e.g., C85; Chakrabarti, 1991) A is given by, 

A 2 = — = + l2gt<t> (20) 

^ 9t<t> + l 9tt 

We choose a flow with equation of state p = KpQ, where K is a constant. Using a distribution 
given by Eq. (19) one obtains the enthalpy h = (p + p)/po = 1/(1 — na 2 s ) (with sound speed a s 
defined by, a 2 s = ^\ s at constant entropy s. n = [7 — 1] _1 is the polytropic index) distribution 
from Eq. (18) as, 

hu t (l - Ql) a = £ , (21) 



« = tt^- (22) 



where, 

rv = 

2n -2 

Equation (21) represents the general relativistic Bernoulli equation. Here, due to the variation of 
angular momentum, the specific energy hu t = £0 is not conserved. Only in an inviscid flow, i.e., 
for n^, = 0, £ = £q. The time component of velocity ut, which is also the specific binding energy, 
could be obtained from Eq. (3) as, 



2 1 1/2 

9t<f> ~ 9tt94>4> 



+ i gt(t) )(i-m)(i-v 2 ) 



(23) 



Before the incoming matter hits this predominantly rotating flow configuration, its thermal 
pressure was negligible, and assuming the same angular momentum distribution (Eq. 19) in the 
pre-shock flow one obtains the ram pressure distribution from the following equation, 

u t (l-ni) a = £, (24) 

Here we have chosen specific enthalpy h = 1 for a cold inflow in the pre-shock region. One need 
not have an angular momentum distribution as given by Eq. (19) in the pre-shock region, since 
the barotropic flow is no longer pre-dominantly rotating, and therefore the von-Zeipel theorem 
(Tassoul, 1978) need not be strictly imposed. But we choose it anyway so as to obtain a purely 
compressible shock wave where the angular momentum remains continuous across the shock. In 
any case, our conclusions are not found to be qualitatively different if the factor (1 — ttl) a is 
dropped instead. 
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At the shock surface, the pressure balance equation (16b) takes the form: 

P+ + (p + + P+ )Vl/(l - Vl) ~ (p_ + P-)^^ (25) 

Note that we have not taken the projection of velocity along the normal to the shock surface here. 
This is because in the strong shock limit that we are considering here, the surfaces of constant ram 
pressure outside the shock and the surfaces of constant thermal pressure inside the shock would 
have similar functional dependence through ut(r, 9). Thus locally the ram pressure gradient is 
parallel to the thermal pressure gradient at each point on the shock surface. 

Combining Eqs. (17a,b) and Eq. (25) we obtain the pressure balance condition at the shock 



as, 



I -rial ~ (7 + l) 2 l-V_ 2 ' [ ' 

Figs, la-b show examples of the shock surfaces (S, solid curves) and the contours (short dashed 
curves) of constant thermal pressure inside the shock surface while the contours of constant ram 
pressure outside the shock surface. Contours are in axial distance R = rsinO vs. vertical distance 
Z = rcosQ (namely, in meridional) plane. The parameters we choose are a = 0.95, = 2.3, 
= 0.05 and 7 = 4/3. We also show the funnel surface defined as the zero pressure surface 
(F) and the centrifugal barrier (CF) defined by the surface of maximum pressure (Vp = 0). In a 
predominantly rotating flow (Eq. 18), the vanishing pressure gradient implies that only competing 
forces at CF are the centrifugal, gravitational and Coriolis' forces (in Newtonian terminology). 
Thus, the centrifugal barrier resembles the Keplerian distribution, so to speak. In Fig. la, we 
choose £ = 1.001 constant at the outer boundary at R = 50. In this case, the shock surface 
that begins on the equatorial plane bends inside. This is because the radial velocity becomes 
smaller as the matter approaches the centrifugal barrier and the flow has to move in farther to 
match the thermal pressure. The segment of the shock surface inside the centrifugal barrier is 
not physical for accretion flow, as the flow has to move outwards to form them. This segment is 
however, physical for winds blowing outwards. Thus both the winds and accretion can have shock 
waves simultaneously in different regions of the flow. The segment outside the centrifugal barrier 
is equivalent to the shock at x s s in the notation of C89 and the segment inside the centrifugal 
barrier is equivalent to x S 2 in the notation of C89. As shown in Chakrabarti & Molteni (1993), in 
a polytropic flow, x s s is stable for accretion while x S 2 is stable for winds. 

In Fig. lb, we artificially choose an outer boundary condition where the specific 
energy of the flow increases as one goes away from the equatorial plane. We choose: 
£{6) = 0.03cos 6 + S{9 = vr/2), where 9 = coincides with the vertical axis. Motivation for 
this choice is that in an astrophysical situation (as well as in a typical numerical simulation) 
infall radial and sound velocities are constants at the outer boundary. At a higher elevation, the 
potential energy decreases, so the net specific energy increases. The flow forms a shock farther 
away from the axis as less centrifugal force is required to balance the weakened gravity. In all 
the numerical simulation results reported (HSW85, HSW85, H84, MLC94, MRC96) such shock 
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waves which bend outwards are seen. Indeed in H84, several examples of shocks which bend back 
towards the vertical axis at a even higher latitude are reported exactly as shown in Fig. lb. After 
the flow passes through this oblique shock, the energy may be randomized, and another shock of 
the type shown in Fig. la may be formed. These are seen in MLC94 and MRC96 as well. 

Note that in the strong shock limit that we discuss here, the outer sonic surface lies at infinity, 
while the inner sonic surface lies arbitrarily close to the horizon. This is simply because the flow is 
chosen to be very cold a s = in the pre-shock flow and almost static V ~ in the post-shock flow. 
When these constraints are relaxed, as in our analysis in the next section, the outer sonic point 
moves in from infinity to a few tens to a few hundreds of Schwarzschild radii, while the inner sonic 
point moves out to roughly one to three Schwarzschild radii as see in pseudo-Newtonian geometry 
(C89, C90). 

4. SHOCKS AND SHOCK-FREE SOLUTIONS IN A CONICAL WEDGE FLOW 

In the last section we concentrated on strong shock solutions. Presently, we study solutions 
of the transonic flows around a rotating black hole in full generality. General procedures for 
obtaining such solutions using pseudo-Newtonian potential have already been in numerous works, 
and we shall not repeat them here (see, C89, C90, C96b). We use this weak viscosity limit of the 
equations of Section 2 to obtain the zeroth order distribution of velocity and sound speed etc. 
Subsequently, we compute the effects of viscosity on this zeroth order distribution. Computation 
of heating and cooling effects are also done in the same way. 

Before proceeding further, we would like to argue that the conical flows have all the properties 
similar to other types of thin flows, such as flows in vertical equilibrium and flows of constant 
thickness. Different models can really be mapped onto each other by varying the polytropic 
constant appropriately. To show this, we consider a gram of matter of specific energy £ and 
specific angular momentum I in Newtonian geometry: 

S= 1 -V 2 + n i a* + ^- 1 -- (27a) 

Using the equation of state p = Kp' 1 (where 7 = 1 + 1/n is the adiabatic index, and K is a 
constant), the mass accretion rate on the equatorial plane 

Mi = pVr q (276) 

could be re- written as, 

Mi = Va?r qi (27c) 

where Mi is an entropy dependent quantity oc MiK ni . Here V is the radial velocity and the 
subscript i differentiates one model from another. For spherical or conical flows, [i = 1, say), 
qi = 2, v\ = 2ni, while for a thin flow of constant height (i = 2, say), 02 = 1, ^2 = 2ri2 and finally, 



- 11 - 



for a thin flow in vertical equilibrium (i = 3, say), 53 = 5/2, ^3 = 2ri3 + 1. It can be easily shown 
that the properties at the sonic points remain identical in all these three models, provided, the 
polytropic indexes are modified according to the relation, 



2ni - 3 2n 2 - 1 2n 3 - 3 



(28) 



Due to non-linearity, a cleaner relation such as this is impossible in general relativity. But we 
expect the situation to be similar. Therefore, we feel confident that a conical wedge shaped flow 
can capture the physical process as any other thin flow model. 

In the case weakly viscous flow, Eq. (6) could be integrated to give the energy conservation, 

P + P 



hu t 



Po 



-u t 



(29a) 



The baryon number conservation relation remains the same as Eq. (7), with the flow cross-section 
area to be A oc 2irr 2 sin 959 where 59 is the angular width of the flow around the angle 9 with 
the vertical axis (in C90, a simpler accretion rate was used to obtain solutions in Kerr geometry) 
Thus, 

M = 2irp r 2 u r = constant (296) 



and finally, the angular momentum conservation equation (Eq. 8) becomes simply, 

P + P 



Po 



(29c) 



We solve these equations by using the standard sonic point analysis. First rewrite the energy and 
mass conservation equations as (e.g., Shapiro & Teukolsky, 1983), 



and 



where, 



and 



M 



I -no 2 (1-^2)1/2 



V 



F(r,9) 



1 — na 2 



(1 - V 2 ) 1 / 2 



F(r,9) 



9t<j> 



G(r,9) 



1/2 



(30a) 



(306) 



r 2 sin9 



G(r,9) = 



1/2 
g rr 



Here, we have neglected 9 component of velocity u$ compared to the other two spatial components 
u r and u^. This assumption need not be valid, especially in the immediate post-shock region, 
where u r itself becomes small enough. However, we believe that we can capture the basic physics 
of transonic flows this approximation. The quantity M will be called the entropy function (earlier, 
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in C89 & C90, the terminology of 'accretion rate' was used for this quantity, though it was 
distinguished from the mass flux M) as Ai oc K n M. It usually assumes two different values at two 
saddle type sonic points and the one with the smallest A4 joins the horizon at infinity (C90). In 
a flow which contains shocks, matter first passes through this transonic solution of lower entropy, 
and then through a shock and finally through the sonic point which requires a higher entropy. 

Differentiating equations (30a) and (30b) with respect to r (for a given 9 = 9q) and eliminating 
terms involving derivatives of a s , we obtain, 

dV N Ci<22 — C2<Zi 



dr D b\a,2 — b 2 a\ 

where, 



(31) 



Zna s 

ai = - ^ ( 32a ) 

I — naj. 

2n 

a 2 = —p. 2V (326) 

a s (i — nag) 

h = (32c) 

d - 1 dFMo) (32e) 

C2 " GMbj dr ' (32/) 

From the vanishing condition of D and N at the sonic points, one obtains the so called sonic 
point conditions as, 

V c = a StC (33a) 

and 

V c 2 = % (336) 

C2 

Here, the subscript 'c' refers to the quantities evaluated at the sonic point r = r c . 

In the study of inviscid accretion flows, one supplies the conserved quantities, such as I, S, 
etc. to obtain a complete solution. For a given angular momentum I, the remaining unknowns are 
V(r), a s (r). But one requires only one extra boundary condition, e.g. £, since two sonic point 
conditions (33a,b) introduce only one extra unknown r c . Both the quantities I and £ could be 
functions of 9 in general and one has to supply another constant equivalent of the Carter constant 
(Chandrasekhar, 1983) to describe flows in the region away from the equatorial plane. However, 
since we let u$ = everywhere, we do not have such a choice any more. Thus, the supply of the 
specific energy 8 and the specific angular momentum I is sufficient to obtain a complete solution 
from the horizon to infinity. Note that by definition Q is the same as uj on the horizon. At 
the shock, apart from the continuity of the energy and the mass flux, one must also satisfy the 
relativistic momentum balance condition, Eq. 16(a-b). 
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4.1. Classification of the Parameter Space 

A solution, with or without a shock, can be completely described by the choice of either (l,M) 
or (l,S) pair. As in a flow in a pseudo-Newtonian geometry (Fig. 2 and Fig. 4 of C89), one can 
classify the entire parameter space in (I, M) or (I, £) plane in terms of whether or not the flow has 
one or more sonic points and shocks. Fig. 2 shows such a classification (for a = 0.95) for the flow 
on the equatorial plane (6 = vr/2). The parameter space is sub-divided into nine regions. Flows 
with parameters from the region N has no transonic solution. The solutions from the region 'O' 
have only the outer sonic point. The solutions from the regions NSA and SA have two 'X' type 
sonic points with the entropy density S at the outer sonic point less than the entropy density Si 
at the inner sonic point. However, flows from SA can pass through a standing shock (See, Figs. 3 
and 4 below) as the Rankine-Hugoniot condition is satisfied. The entropy generated at the shock 
is exactly Si — S Q , which is advected towards the black hole to enable the flow to pass through 
the inner sonic point. Rankine-Hugoniot condition is not satisfied for flows from the region NSA. 
Numerical simulation indicates (Ryu, Chakrabarti & Molteni, 1997) that the solution is very 
unstable and show periodic changes in shock locations. The flow from the region SW and NSW 
are very similar to those from SA and NSA. However, S > Si in these cases. Shocks can form 
only in winds from the region SW. The shock condition is not satisfied in winds from the region 
NSW. This may make the NSW flows unstable as well. However, stable accretion can take place 
through the inner sonic point. A flow from region / only has the inner sonic point and thus can 
form shocks (which require the presence of two saddle type sonic points) only if the inflow is 
already supersonic due to some other physical processes. See C90 and Chakrabarti (1996c) for all 
the topologies of the flow. 

The flows from regions I* and O* are interesting in the sense that each of them has two sonic 
points (one 'X' type and one 'O' type) only and neither produces a globally complete solution 
(in the inviscid case, at least). The region I* has an inner sonic point but the solution does not 
extend subsonically to a large distance. The region O* has an outer sonic point, but the solution 
does not extend supersonically to the horizon! In both the cases a weakly viscous flow is expected 
to be unstable. When a significant viscosity is added, the closed topology of /* is expected to 
open up (C90, C96b) and then the flow can join with a Keplerian disk. Clearly, advective disks 
from this region of parameter space will have no standing shocks independent of what polytropic 
index is used unless some extra energizing mechanism (such as flares) in the sub-Keplerian region 
raises the specific energy of the flow to bring it to the parameter space marked 'SA'. 

4.2. Behavior of Shock Locations on the Equatorial Plane 

First, we study the properties of the shock waves on the equatorial plane {9 = n/2). Our 
goal is to see if the shocks in a true black hole geometry display similar properties as those in 
pseudo-Newtonian geometry. With a pseudo-Newtonian geometry, it was noted that for given 
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flow parameters, there are two locations x S 2 and x s z where shocks form in between two saddle 
type sonic points. Later, it was found (Chakrabarti & Molteni, 1993; also see Nakayama, 1992 for 
isothermal flows) that x s ^ is stable in accretion and x S 2 is stable in winds. In Fig. 2, we show 
shock locations in four cases, (a) a = 0, I = 3.5, (b) a = 0.5, I = 3.0, (c) a = 0.95, I = 2.3 and (d) 
a = 0.95, I = 2.1 as a function of the specific energy of the flow. Each 'Persian bottle' shaped set 
is marked with the above parameters. Each set contains four curves, two dashed and two solid. 
The dashed curve on the lower left (02) corresponds to the unstable shock x S 2 of accretion flow 
and the solid curve on the upper left (03) corresponds to the stable shock x S 3 of accretion flow. 
The solid curve on the lower right (^2) corresponds to the stable shock x S 2 of the wind flow and 
the dashed curve on the upper right (103) corresponds to the unstable shock x s s of wind flow. This 
description is applicable to each set. First notice that the shocks are possible on the equatorial 
plane for a certain range of energy values only. The general behavior of the shock location is 
exactly same as in pseudo- Newtonian geometry (C89; C90; Chakrabarti & Molteni, 1993). What 
is more, the range of energy for which shocks are possible in a rapidly rotating black hole is much 
higher compared to that around a non-rotating black hole. On the average, the shock locations 
are closer to the horizon for flows around a rotating black hole. 

4.3. Shear, Heating and Cooling in a Viscous flow 

Since the nature of the transonic solution including the formation of shocks depends on 
angular momentum distribution, which in turn depend very strongly upon viscosity cooling and 
heating, it is instructive to study their behavior. Presently, we study then in the weak viscosity 
limit, where the velocity and density distributions are chosen from the inviscid flows while 
computing shear, heating and cooling. The complete solutions of the viscous transonic flow which 
self-consistently compute these quantities along with the flow variables (such as, velocity, density, 
etc.) will be dealt with elsewhere. 

Around a Newtonian star, shear component due to purely rotational motion is defined 
as a r( p = rdfl/dr which is always positive (e.g, 3/20 for a Keplerian disk), and hence viscous 
stress is always negative. That is, angular momentum is always transported outwards. However, 
Anderson h Lemos (1988) pointed out that the shear can change sign near the horizon and they 
demonstrated this using velocity profiles of a cold radial flow below marginally stable orbit. Below, 
we show that at least for a weakly viscous flow shear indeed changes its sign close to a black hole. 

In Fig. 4, we show the behavior of shear components crj for a prograde and a retrograde 
flow. In the upper left panel, we show the Mach number variation of the flow with 
a = 0.95, I = 2.3, £ = 1.001. The arrowed path shows that two transonic solutions are 'glued' 
together by the vertical shock transition (at 03 or x S 3 in our notation) where Rankine Hugoniot 
condition is satisfied. On the upper right panel, We present a'^ p and cr'^\ ro t which are the 
general shear and the rotational shear (Eq. 9) for the supersonic branch passing through the 
outer sonic point (subscripts 'p' and 'b' are for supersonic and subsonic branches, respectively). 
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For comparison, we also show a r , h (general shear for the subsonic branch passing through the 
outer sonic point) and dVi/dr (which is independent of any branch). One arrives at several 
important conclusions: First, the effect of radial motion is significant in changing the shear stress. 
The rotational shear <T^| rot as well as a^ b where the radial velocity is negligible, vanish at the 
horizon, while the a^ p does not. Indeed, a'^ p is a very large negative quantity, and therefore in 
presence of even a weak viscosity it can transport angular momentum inwards. This does not 
automatically imply that the angular momentum distribution will develop a minimum (such as in 
the Keplerian distribution), see Fig. 5 below. In reality, this inward transport may try to nullify 
a steeper distribution with a positive slope and makes it as flat as possible near the horizon. 
The lower left panel of Fig. 4 gives the variation of Mach numbers for a retrograde flow with 
a = 0.95, I = —4.0, £ = 1.0045. Note that both the subsonic and supersonic branches behave 
qualitatively similarly, although that from the supersonic branch (cr r ^ p ) becomes extremely large 
positive quantity before turning back to become a large negative just outside the horizon. We have 
not plotted the shear distribution for the post-shock branches passing through the inner sonic 
points, as their behaviors are similar. What is very much clear is that a simple Shakura-Sunyaev 
(1973) type viscous stress oc — ap with a constant a will not be very accurate, specially close to 
the black hole, since the pressure p is always positive while actual stress changes sign! 

The reversal of the shear close to a black hole is no mystery. In general relativity, all the 
energies couple one another. It is well known that the 'pit in the potential' of a black hole is due 
to coupling between the rotational energy and and the gravitational energy (see, Chakrabarti, 
1993 and references therein). As matter approaches the black hole, the rotational energy, and 
therefore 'mass' due to the energy increases which is also attracted by the black hole. This makes 
gravity much stronger than that of a Newtonian star. When the black hole rotates, there are 
more coupling terms (such as that arises out of spin of the black hole and the orbital angular 
momentum of the matter) which either favor gravity or go against it depending on whether the 
flow is retrograde or prograde respectively. This is the basic reason why retrograde and prograde 
flows display different reversal properties. 

Given that the shear is a very large negative number close to a horizon, for a physical flow, 
one must have the viscosity vanishing on the horizon in order to have finite angular momentum 
or the viscous dissipation. Fig. 5 shows the distribution of angular momentum (upper panel) and 
the heating and cooling rates (lower panel) in the pre-shock (with subscript 'pre') and post-shock 
(with subscript 'post') solutions for the prograde flow presented in Fig. 4 above. We also present 
the Keplerian distribution (Ik) in the upper panel for comparison. We assume ion viscosity 
(Spitzer, 1962) with a factor (1 — r + /r) lS in order that it vanishes on the horizon sufficiently 
rapidly (we use (3 = 1 here). The same variation is used to compute viscous heating rate: 

Q + = 2i 1 o^ v = -^a r(/) al(n 2 g^ + 2Slg t(j> + g tt ) (34) 

Note that the quantity in the parenthesis has to be negative in order to obtain positive heating. 
This is indeed always the case for any realistic flow (e.g., Shapiro & Teukolsky, 1983). For the 
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cooling, we simply use bremsstrahlung process. In presence of a Keplerian disk farther out, 
Comptonization of the intercepted photons is likely to dominate (Chakrabarti & Titarchuk, 
1995). We chose = 1M and M = MEddington f° r illustration purposes (as we are interested 
in general behavior, and not the absolute magnitude). In the upper panel, dotted curve is the 
angular momentum distribution in the preshock branch and solid curve is that in the post-shock 
branch. The dotted vertical arrow is at the shock location. In the lower panel the solid curves 
are viscous dissipation rate and the dotted curves are the bremsstrahlung cooling rate. Clearly, 
since the post-shock flow is hotter, ion viscosity is higher. This causes a rapid transfer of angular 
momentum as well as a higher heating and cooling rates in the post-shock region. In reality, when 
flow velocity, sound speed etc. are computed simultaneously, the effect becomes much weaker, 
and the piling up of angular momentum would either move the shock outward, or the shock 
may disappear completely as has been observed in recent numerical simulations (see, C96a and 
References therein). Most of the heat generated in the present case is carried away towards the 
black hole, since bremsstrahlung is not efficient enough to cool the flow. This is similar to the 
inviscid case where energy remained conserved (C89), or in viscous pseudo-Newtonian flows with 
inefficient cooling (C96a, C96b). 

Instead of using ion viscosity, if one used parametrized kinematic viscosity such as v ~ ac s r 
where, c s is the isothermal sound speed and a is the Shakura-Sunyave (1973) viscosity parameter, 
then also we would have exactly the same behavior as above. However, for finite heating and 
angular momentum transport rate, a must vanish 'sufficiently rapidly' on the horizon. This 
suggests that it is essential to have a fundamental change in viscous processes very close to the 
horizon. A constant a description may be adequate for 'predominately rotating flows', however, 
since the rotational shear <r r ,\ ro t vanishes on the horizon (Fig. 4). Such predominantly rotating 
flows do not satisfy the boundary condition V = 1 on the horizon and we prefer to conjecture that 
the turbulence and viscosity be completely damped out on the horizon instead ('tranquil inflow 
conjecture'). 



4.4. Shock Location away from the Equatorial plane 

If one considers conical flows away from the equatorial plane, the shocks are generally 
expected to recede from the axis. This is simply because away from the plane the gravitational 
force diminishes so that less centrifugal force is required to form a shock. For illustration, we 
consider the same three cases as in the above sub-section, but we choose a particular energy 
£ = 1.004 for all the three cases. From Fig. 3, we note that shock solutions on the equatorial 
plane exist for all the three cases at this energy. In Fig. 6, we show the variation of shock location 
(dotted curves marked with S) with height in the R — Z plane. For comparison, we show the 
funnel wall (solid curves marked with F) and the centrifugal barrier (long dashed curves marked 
with CF). The outer most curve in each set is for a = 0, the middle curve is for a = 0.5 and the 
innermost curve is for a = 0.95. The location of the shock is more distant from the vertical axis at 
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a higher elevation. The shock strengths (measured by the ratio of pre-shock to post-shock Mach 
numbers M_/M + ) fall sharply with increasing latitude. The shocks in accretion flows typically 
exist in the range 9 max > 9 > 9 m i na (see Fig. 7 below), where 9 max is not necessarily ir/2. In 
other words, even for parameters which do not allow shocks on the equatorial plane, the flow may 
have shocks starting at an angle 9 max < vr/2 away from the equatorial plane and extending up to 
9min,a- This happens when the energy is less than the lower limit of energy required for the shock 
in the equatorial plane (which is different for different cases as depicted in Fig. 2 and Fig. 3). For 
instance, flows from the region NSA (Fig. 2) would have shocks at a higher elevation, though on 
the equatorial plane shocks were not allowed. Beyond 9 m i nA , the flow can have stable shock x S 2 
only in winds up to some other 9 m i n ^ w . These shocks form very close to the centrifugal barrier 
near the black hole, and for clarity we do not draw them. These shocks also bend away from the 
vertical axis at higher elevations. Similar to accretion flows, winds can develop shocks at higher 
elevation when the energy is more than the highest energy required for shocks on the equatorial 
plane (which is different in different cases as depicted in Fig. 2 and Fig. 3). Thus flows from 
region NSW (Fig. 2) would have shocks in winds at higher elevation. Again for clarity we have 
not plotted outer and inner sonic surfaces. These 'annular shocks' are observed in simulations of 
Ryu et al. (1995) and Ryu, Chakrabarti and Molteni, 1996). The outer sonic surface is very close 
to spherical (as it is hundreds of Schwarzschild radii away, and therefore insensitive to rotational 
effects) moving in slightly with elevation. The inner sonic surface located close to the horizon 
moves inwards with elevation at a comparatively faster rate. 

In Fig. 7 we compare the entropy function M for all the three cases discussed in §4.2 as a 
function of 9 (9 = 90 corresponds to the equatorial plane). For a given mass accretion rate M for 
all the three cases, these curves measure the entropy of the flow and the amount of entropy that is 
generated at the shock. The three dashed curves, roughly overlapping each other are the entropy 
functions at the outer sonic surfaces. These sonic surfaces being very far away, their properties are 
almost independent of the state of the black hole or the (sub-Keplerian) angular momentum of 
matter. These are the entropies of the pre-shock flow since no entropy is generated before the flow 
hits the shock. The solid curves represent the entropy of the post-shock flow which pass through 
the inner sonic surfaces. Clearly, the flow around the rapidly rotating black hole is the hottest and 
has the highest entropy. The difference in entropy between the outer and the inner sonic surfaces 
must be generated at the shock. Therefore, a higher jump (at a given angle) in entropy represents 
a stronger shock. Note that at some 9 m i n ^ a the solid and dashed curves intersect. That is the 
location (angle) where the weakest shock (with strength unity) forms. For 9 m i n ^ a < 9 < ir/2, the 
shock in accretion flow can form, since it has to generate entropy at the shock. Beyond that the 
entropy at the outer sonic surface is lower than that from the inner sonic surface, and the shocks 
are formed only in winds. Accretion flow with these parameters will not have shocks. This process 
continues till # m j n ,«> ~ 60 degrees. Beyond that angle, shocks will not form in winds as well. 
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5. CONCLUDING REMARKS 

In this paper, we studied the general behavior of the transonic solutions, with or without 
shocks, in two dimensional accretion and winds. We employed fully general relativistic equations 
for viscous transonic flows for this purpose. We show that the behavior of shock waves on the 
equatorial plane is very similar to what had been found using pseudo-Newtonian potential. 
However, for rotating black holes, the range of parameter space which allows the formation of 
shock waves is much larger. Similarly, flows away from the plane form shocks more easily. These 
studies therefore imply that shocks in quasi-spherical accretion flows are perhaps more generic. 
Furthermore, it can be easily shown, by employing a perturbation of pressures on both sides of 
the shocks, that these shocks are stable as well (Chakrabarti & Molteni, 1993). Both of these 
deductions from purely analytical considerations as presented in this paper are borne out by 
every numerical simulations to date (HSW84, HSW85, H84, MLC94, MRC96). These simulations 
always find very large scale stable oblique shocks in a black hole accretion. In case when the flow 
is bound, namely £ < 1, the flow does not have outer sonic point, and does not extend to a large 
distance (region I* in Fig. 2). Only when viscosity is added, the closed topology around 'O' type 
sonic point opens up to join a Keplerian disk farther out. Shocks can form on the equatorial plane 
in this case only when the inflow is already supersonic, though away from the plane the flow may 
have £ > 1 and stable shocks will form. This considerations were used to explain the spectral 
properties of black hole candidates (Chakrabarti & Titarchuk, 1995). In the advective corona 
£ > 1 condition could also be achieved by energy deposition in the corona through magnetic flares 
(for instance) a method which is widely accepted in high energy astrophysics. 

Our 'brute force' method of global shock study (Section 3) solves the mystery that was first 
noticed by Fukue (1987) and subsequently by C89 as to the multiplicity in formal shock locations. 
Since the transonic flow we consider are necessarily non-Kepler ian, they typically tend to have 
a pressure maximum exactly as in a thick accretion disk. As a result, the shock condition is 
satisfied once on each side of the maximum. Our method also reveals why the outer shock is 
stable in accretion and inner shock is stable in winds. In both the cases we see that the thermal 
pressure must increase downstream. In fact, for a highly viscous flow the pressure maxima and 
consequently the shocks are absent (C96b). 

Since general topology of the solutions depend on angular momentum distribution, we also 
studied the behavior of shear for a transonic flow. We conclude that shear is not monotonic close 
to a black hole. In fact it may reverse several times due to coupling of various energies in general 
relativity. Because of this we believe that neither a viscosity prescription, nor computation based 
on predominantly rotating flow (rotational shear, Eq. [9]) should be adequate. We find that 
heating is very strong close to the black hole, and the generated heat must be advected away with 
the flow if the cooling is not very efficient. This is similar to inviscid case (e.g. C89), where the 
entire energy and entropy at the shock was advected by the optically thin flow. The post-shock 
solution is observed to transport angular momentum more rapidly and as a result the shock may 
either move away from the black hole, or it may disappear if the viscosity is high enough. These 
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results agree with our earlier computations in pseudo-Newtonian flows. 

Since black holes have no hard surface one would imagine that matter would enter through 
the horizon without decelerating. But it seems that the centrifugal barrier by even a small amount 
of angular momentum in the flow causes the matter to pile up behind the barrier over a large 
scale, not just near the equatorial plane. Indeed, we found that the probability of formation of 
the shocks away from the plane is much higher, simply because the barrier is definitely present 
away from the plane (due to weaker gravity) even when it is absent on the equatorial plane (such 
as when the angular momentum of the flow is everywhere sub-Keplerian close to the black hole). 
Our present findings support the view that the black hole accretion models may need to include 
shock waves (at least for flows at higher elevation) as they provide a complete explanation of 
the observed steady state as well as time dependent behaviors (e.g., Chakrabarti & Titarchuk, 
1995; Ryu, Chakrabarti & Molteni, 1997; Crary et al, 1996). For instance, according to our 
model, quasi-periodic oscillations of black hole candidates should have frequencies of the order of 
l/(4x^ 2 ) (Molteni, Sponholz & Chakrabarti, 1996) and therefore can vary from a few millihertz 
to kilohertz range depending on shock locations (Fig. 3) with modulation amplitude as much as 
100 percent. This roughly agrees with the observed frequencies and amplitudes. 

The author is thankful to his colleagues D. Ryu and D. Molteni for excellent numerical 
simulations which partly motivated him to complete the present work. He is also thankful to J. 
Hawley for giving him a copy of his thesis ten years ago, which contained insightful illustrations 
which were helpful in understanding two dimensional transonic flows. He thanks J. Peitz, R. 
Khanna and Kip Thorne for discussions. 
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FIGURE CAPTIONS 

Fig. la-b Funnel wall (F), centrifugal barrier (CF) and the strong shock surface (S) in a flow 
which is cold in the pre-shock region and predominantly rotating in the post-shock region. 
The dashed contours inside the shock are those of thermal pressure plus ram pressure while 
those outside are those of ram pressure. In regions outside the centrifugal barrier, the shock 
is stable in accretion while that inside is stable in winds. In (a) energy is constant at a 
constant von-Zeipel cylindrical surface, while in (b) energy is increasing with latitude to 
mimic boundary conditions of numerical simulations. 

Fig. 2 Classification of the entire parameter space spanned by the specific energy and angular 
momentum in terms of the number of sonic points and the presence or absence of shocks. 
See text for details. 

Fig. 3 Stable (solid) and unstable (dashed) shock locations in accretion (tt2, 03) and winds 
(w2 and W3) as a function of the specific energy energy of the flow for four black holes: 
a = 0, 1 = 3.5; a = 0.5, I = 3.0; a = 0.95, I = 2.3; a = 0.99, I = 2.1. 

Fig. 4 Variation of Mach numbers (left panels) and shear components (right panels) in a prograde 
(upper panels) and retrograde (lower panels) flow. The shock transition is shown as a vertical 
arrow. Shear components are not monotonic as in a Newtonian flow, but reverses near the 
black hole. o r ^\rot is for rotational shear (assuming the flow is predominantly rotating), 
whereas a r , ^ b are for general shear which includes the radial motion. 

Fig. 5 Variation of the angular momentum distribution (upper panel) and the heating and cooling 
rates (lower panel) in presence of ion- viscosity and bremsstrahlung effects respectively. Most 
of the heat generated by the viscous transonic flow in this case is carried by the flow inwards 
as the bremsstrahlung is inefficient in radiating it away. 

Fig. 6 Funnel walls (Solid curves, marked by F), centrifugal barriers (long dashed curves, marked 
by CF) and shock surfaces (short dashed curves, marked by S) for the three cases of Fig. 
2 for £ = 1.004. In each set, the outermost, middle and the innermost curves are for 
a = 0, 0.5, 0.95 respectively. 

Fig. 7 Entropy function Ad is plotted against the inclination angle 9 for the three cases as in Fig. 
3. The dashed curved are the values for flows passing through the outer sonic surfaces A4 , 
and the solid curves are for flows passing through the the inner sonic surfaces M{. Shocks 
in accretion are possible only when A4i > A4 a , namely when Q m in,a < < 6 m ax = vr/2 
(equatorial plane). m in,a is ~ 79°, ~ 72° and 70° respectively for a = 0, 0.5, 0.95 
respectively. For ~ 60°^# < 9 m in,a shocks may form in winds. 
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